# This file produces figures in the text and supplementary information for
#   Goehring and Lowande's Public Responses to Unilateral Policymaking. 

# See the README for more information. 
require(arm)
require(lubridate)
require(gridExtra)
require(ordinal)
require(stargazer)
require(kableExtra)
require(mfx)
require(xtable)
require(slider)
require(tidyverse)
require(broom)
require(irr)
require(MASS)

load('input-data/complete-data-v3.Rda') 


# TOPIC APPROVAL
t1_d=lm(d_topic~topic.w2+president.w2+condition.w2+outcome,data=dt)
t2_d=lm(d_topic~topic.w2+president.w2+condition.w2+partisanship+sex+age+income+ethnicity+outcome,data=dt)
t3_d=lm(d_topic~topic.w2+president.w2+condition.w2+partisanship+sex+age+income+ethnicity+education+outcome,data=dt)

# PRESIDENTIAL APPROVAL
p1_d=lm(d_job~topic.w2+president.w2+condition.w2+outcome,data=dt)
p2_d=lm(d_job~topic.w2+president.w2+condition.w2+partisanship+sex+age+income+ethnicity+outcome,data=dt)
p3_d=lm(d_job~topic.w2+president.w2+condition.w2+partisanship+sex+age+income+ethnicity+education+outcome,data=dt)

# TRUMP VOTE
v1_d=lm(d_trump~topic.w2+president.w2+condition.w2+outcome,data=dt)
v2_d=lm(d_trump~topic.w2+president.w2+condition.w2+partisanship+sex+age+income+ethnicity+outcome,data=dt)
v3_d=lm(d_trump~topic.w2+president.w2+condition.w2+partisanship+sex+age+income+ethnicity+education+outcome,data=dt)

# GET DONE
g1_d=lm(d_done~topic.w2+president.w2+condition.w2+outcome,data=dt)
g2_d=lm(d_done~topic.w2+president.w2+condition.w2+partisanship+sex+age+income+ethnicity+outcome,data=dt)
g3_d=lm(d_done~topic.w2+president.w2+condition.w2+partisanship+sex+age+income+ethnicity+education+outcome,data=dt)

# condition on pre-treatment value

# TOPIC APPROVAL
t1_c=clm(topic.approve.w2~topic.w2+president.w2+condition.w2+topic.approve.w1+outcome,data=dt,family='logit')
t2_c=clm(topic.approve.w2~topic.w2+president.w2+condition.w2+topic.approve.w1+partisanship+sex+age+income+ethnicity+outcome,data=dt,family='logit')
t3_c=clm(topic.approve.w2~topic.w2+president.w2+condition.w2+topic.approve.w1+partisanship+sex+age+income+ethnicity+education+outcome,data=dt,family='logit')

# PRESIDENTIAL APPROVAL
p1_c=clm(pres.approve.w2~topic.w2+president.w2+condition.w2+pres.approve.w1+outcome,data=dt,family='logit')
p2_c=clm(pres.approve.w2~topic.w2+president.w2+condition.w2+pres.approve.w1+partisanship+sex+age+income+ethnicity+outcome,data=dt,family='logit')
p3_c=clm(pres.approve.w2~topic.w2+president.w2+condition.w2+pres.approve.w1+partisanship+sex+age+income+ethnicity+education+outcome,data=dt,family='logit')

# TRUMP VOTE
v1_c=glm(trumpvote.w2~topic.w2+president.w2+condition.w2+trumpvote.w1+outcome,data=dt,family='binomial')
v2_c=glm(trumpvote.w2~topic.w2+president.w2+condition.w2+trumpvote.w1+partisanship+sex+age+income+ethnicity+outcome,data=dt,family='binomial')
v3_c=glm(trumpvote.w2~topic.w2+president.w2+condition.w2+trumpvote.w1+partisanship+sex+age+income+ethnicity+education+outcome,data=dt,family='binomial')

# GET DONE
g1_c=clm(manip.done.w2~topic.w2+president.w2+condition.w2+manip.done.w1+outcome,data=dt,family='logit')
g2_c=clm(manip.done.w2~topic.w2+president.w2+condition.w2+manip.done.w1+partisanship+sex+age+income+ethnicity+outcome,data=dt,family='logit')
g3_c=clm(manip.done.w2~topic.w2+president.w2+condition.w2+manip.done.w1+partisanship+sex+age+income+ethnicity+education+outcome,data=dt,family='logit')

# rerun models with binary DV
t3_c.bin=glm(topic.approve.w2.bin~topic.w2+president.w2+condition.w2+topic.approve.w1+partisanship+sex+age+income+ethnicity+education+outcome,data=dt,family='binomial')
p3_c.bin=glm(pres.approve.w2.bin~topic.w2+president.w2+condition.w2+pres.approve.w1+partisanship+sex+age+income+ethnicity+education+outcome,data=dt,family='binomial')
g3_c.bin=glm(manip.done.w2.bin~topic.w2+president.w2+condition.w2+manip.done.w1+partisanship+sex+age+income+ethnicity+education+education+outcome,data=dt,family='binomial')

# function for simulating marginal changes
outcome_comparisons=function(model,data,simulations=1000,seed=6658) {
  set.seed(seed)
  betas=slot(sim(model,simulations),'coef')
  covariates=names(model$model)
  results=matrix(data=NA, ncol=7,nrow=1)
  
  # add columns for factor coefs
  D=data[complete.cases(data[,covariates]),]
  D=cbind(D,model.matrix(~topic.w2-1,data=D))
  D=cbind(D,model.matrix(~president.w2-1,data=D))
  if ('topic.approve.w1' %in% covariates) {D=cbind(D,model.matrix(~topic.approve.w1-1,data=D))}
  if ('pres.approve.w1' %in% covariates) {D=cbind(D,model.matrix(~pres.approve.w1-1,data=D))}
  if ('manip.done.w1' %in% covariates) {D=cbind(D,model.matrix(~manip.done.w1-1,data=D))}
  if ('partisanship' %in% covariates) {D=cbind(D,model.matrix(~partisanship-1,data=D))}
  if ('sex' %in% covariates) {D=cbind(D,model.matrix(~sex-1,data=D))}
  if ('income' %in% covariates) {D=cbind(D,model.matrix(~income-1,data=D))}
  if ('ethnicity' %in% covariates) {D=cbind(D,model.matrix(~ethnicity-1,data=D))}
  if ('education' %in% covariates) {D=cbind(D,model.matrix(~education-1,data=D))}
  D=cbind(D,model.matrix(~condition.w2-1,data=D))
  
  # intercept + covariates 
  controls=cbind(1,D$`topic.w2contractor pay`,D$`topic.w2farm payments`,D$`topic.w2foreign aid and abortion`,D$`topic.w2foreign workers`,D$`topic.w2gun violence research`,D$topic.w2LGBT,D$`topic.w2policing and military surplus`,D$`topic.w2public lands`,D$`topic.w2Russian sanctions`,D$`topic.w2student loans`,D$topic.w2trade,D$`topic.w2water rules`,D$`topic.w2weapon sales to Saudi Arabia`,D$topic.w2wildlife,D$president.w2Trump,D$condition.w2congress,D$condition.w2order)
  if ('topic.approve.w1' %in% covariates) {controls=cbind(controls,D$topic.approve.w1Disapprove,D$`topic.approve.w1Somewhat disapprove`,D$`topic.approve.w1Neither approve or disapprove`,D$`topic.approve.w1Somewhat approve`,D$`topic.approve.w1Approve`,D$`topic.approve.w1Strongly approve`)}
  if ('pres.approve.w1' %in% covariates) {controls=cbind(controls,D$pres.approve.w1Disapprove,D$`pres.approve.w1Somewhat disapprove`,D$`pres.approve.w1Neither approve or disapprove`,D$`pres.approve.w1Somewhat approve`,D$`pres.approve.w1Approve`,D$`pres.approve.w1Strongly approve`)}
  if ('manip.done.w1' %in% covariates) {controls=cbind(controls,D$manip.done.w1Disagree,D$`manip.done.w1Somewhat disagree`,D$`manip.done.w1Neither agree nor disagree`,D$`manip.done.w1Somewhat agree`,D$`manip.done.w1Agree`,D$`manip.done.w1Strongly agree`)}
  if ('trumpvote.w1' %in% covariates) {controls=cbind(controls,D$trumpvote.w1)}
  if ('partisanship' %in% covariates) {controls=cbind(controls,D$partisanshipDemocrat,D$partisanshipRepublican)}
  if ('sex' %in% covariates) {controls=cbind(controls,D$sexFemale)}
  if ('age' %in% covariates) {controls=cbind(controls,D$age)}
  if ('income' %in% covariates) {controls=cbind(controls,D$`income2: $25,000 to $50,000`,D$`income3: $50,001 to $75,000`,D$`income4: $75,001 to $100,000`,D$`income5: More than $100,001`)}
  if ('ethnicity' %in% covariates) {controls=cbind(controls,D$ethnicityBlack,D$`ethnicityNative American`,D$`ethnicityOther/Decline to State`,D$ethnicityWhite)}
  if ('education' %in% covariates) {controls=cbind(controls,D$`educationSome College or Vocational`,D$`educationB.A. or B.S.`,D$`educationPost-graduate or Higher`)}
  
  # new data to simulate over
  d.success=cbind(controls,0)
  d.fail=cbind(controls,1)
  
  # simulations
  p.win=apply((1/(1+exp(-as.matrix(d.success)%*% t(betas)))),1,as.vector) 
  p_win=apply(p.win,1,mean)
  p.fail=apply((1/(1+exp(-as.matrix(d.fail)%*% t(betas)))),1,as.vector) 
  p_fail=apply(p.fail,1,mean)
  
  # store
  results[1,1]<- quantile(p_fail-p_win,.0083)
  results[1,2]<- quantile(p_fail-p_win,.025)
  results[1,3]<- mean(p_fail-p_win)
  results[1,4]<- quantile(p_fail-p_win,.975)
  results[1,5]<- quantile(p_fail-p_win,.9916)
  results[1,6]<- 'All Respondents'
  results[1,7]<- covariates[1]
  
  return(results)
}

outcomes=data.frame(lb.MC=numeric(4),lb.95=numeric(4),me=numeric(4),ub.95=numeric(4),ub.MC=numeric(4),grp=character(4),dv=character(4))
outcomes[1,1:7]=outcome_comparisons(t3_c.bin,dt)
outcomes[2,1:7]=outcome_comparisons(p3_c.bin,dt)
outcomes[3,1:7]=outcome_comparisons(v3_c,dt)
outcomes[4,1:7]=outcome_comparisons(g3_c.bin,dt)

# copartisanship interactions

# change score

# TOPIC APPROVAL
t1a_d=lm(d_topic~topic.w2+president.w2+condition.w2+copartisanship*outcome,data=dt)
t2a_d=lm(d_topic~topic.w2+president.w2+condition.w2+sex+age+income+ethnicity+copartisanship*outcome,data=dt)
t3a_d=lm(d_topic~topic.w2+president.w2+condition.w2+sex+age+income+ethnicity+education+copartisanship*outcome,data=dt)

# JOB APPROVAL
p1a_d=lm(d_job~topic.w2+president.w2+condition.w2+copartisanship*outcome,data=dt)
p2a_d=lm(d_job~topic.w2+president.w2+condition.w2+sex+age+income+ethnicity+copartisanship*outcome,data=dt)
p3a_d=lm(d_job~topic.w2+president.w2+condition.w2+sex+age+income+ethnicity+education+copartisanship*outcome,data=dt)

# TRUMP VOTE
v1a_d=lm(d_trump~topic.w2+president.w2+condition.w2+copartisanship*outcome,data=dt)
v2a_d=lm(d_trump~topic.w2+president.w2+condition.w2+sex+age+income+ethnicity+copartisanship*outcome,data=dt)
v3a_d=lm(d_trump~topic.w2+president.w2+condition.w2+sex+age+income+ethnicity+education+copartisanship*outcome,data=dt)

# GET DONE
g1_d=lm(d_done~topic.w2+president.w2+condition.w2+copartisanship*outcome,data=dt)
g2_d=lm(d_done~topic.w2+president.w2+condition.w2+sex+age+income+ethnicity+copartisanship*outcome,data=dt)
g3_d=lm(d_done~topic.w2+president.w2+condition.w2+sex+age+income+ethnicity+education+copartisanship*outcome,data=dt)

# condition on pre-treatment value

# TOPIC APPROVAL
t1a_c=clm(topic.approve.w2~topic.w2+president.w2+condition.w2+topic.approve.w1+copartisanship*outcome,data=dt,family='logit')
t2a_c=clm(topic.approve.w2~topic.w2+president.w2+condition.w2+topic.approve.w1+sex+age+income+ethnicity+copartisanship*outcome,data=dt,family='logit')
t3a_c=clm(topic.approve.w2~topic.w2+president.w2+condition.w2+topic.approve.w1+sex+age+income+ethnicity+education+copartisanship*outcome,data=dt,family='logit')

# JOB APPROVAL
p1a_c=clm(pres.approve.w2~topic.w2+president.w2+condition.w2+pres.approve.w1+copartisanship*outcome,data=dt,family='logit')
p2a_c=clm(pres.approve.w2~topic.w2+president.w2+condition.w2+pres.approve.w1+sex+age+income+ethnicity+copartisanship*outcome,data=dt,family='logit')
p3a_c=clm(pres.approve.w2~topic.w2+president.w2+condition.w2+pres.approve.w1+sex+age+income+ethnicity+education+copartisanship*outcome,data=dt,family='logit')

# TRUMP VOTE
v1a_c=glm(trumpvote.w2~topic.w2+president.w2+condition.w2+trumpvote.w1+copartisanship*outcome,data=dt,family='binomial')
v2a_c=glm(trumpvote.w2~topic.w2+president.w2+condition.w2+trumpvote.w1+sex+age+income+ethnicity+copartisanship*outcome,data=dt,family='binomial')
v3a_c=glm(trumpvote.w2~topic.w2+president.w2+condition.w2+trumpvote.w1+sex+age+income+ethnicity+education+copartisanship*outcome,data=dt,family='binomial')

# GET DONE
g1a_c=clm(manip.done.w2~topic.w2+president.w2+condition.w2+manip.done.w1+copartisanship*outcome,data=dt,family='logit')
g2a_c=clm(manip.done.w2~topic.w2+president.w2+condition.w2+manip.done.w1+sex+age+income+ethnicity+copartisanship*outcome,data=dt,family='logit')
g3a_c=clm(manip.done.w2~topic.w2+president.w2+condition.w2+manip.done.w1+sex+age+income+ethnicity+education+copartisanship*outcome,data=dt,family='logit')

# rerun with binary DV for simulations
t3a_c.bin=glm(topic.approve.w2.bin~topic.w2+president.w2+condition.w2+topic.approve.w1+sex+age+income+ethnicity+copartisanship*outcome,data=dt,family='binomial')
p3a_c.bin=glm(pres.approve.w2.bin~topic.w2+president.w2+condition.w2+pres.approve.w1+sex+age+income+ethnicity+copartisanship*outcome,data=dt,family='binomial')
g3a_c.bin=glm(manip.done.w2.bin~topic.w2+president.w2+condition.w2+manip.done.w1+sex+age+income+ethnicity+copartisanship*outcome,data=dt,family='binomial')

outcome_comparisons_party=function(model,data,simulations=1000,seed=6658) {
  set.seed(seed)
  betas=slot(sim(model,simulations),'coef')
  covariates=names(model$model)
  results=matrix(data=NA, ncol=7,nrow=3)
  
  # add columns for factor coefs
  D=data[complete.cases(data[,covariates]),]
  D=cbind(D,model.matrix(~topic.w2-1,data=D))
  D=cbind(D,model.matrix(~president.w2-1,data=D))
  if ('copartisanship' %in% covariates) {D=cbind(D,model.matrix(~copartisanship-1,data=D))}
  if ('topic.approve.w1' %in% covariates) {D=cbind(D,model.matrix(~topic.approve.w1-1,data=D))}
  if ('pres.approve.w1' %in% covariates) {D=cbind(D,model.matrix(~pres.approve.w1-1,data=D))}
  if ('manip.done.w1' %in% covariates) {D=cbind(D,model.matrix(~manip.done.w1-1,data=D))}
  if ('copartisanship' %in% covariates) {D=cbind(D,model.matrix(~partisanship-1,data=D))}
  if ('sex' %in% covariates) {D=cbind(D,model.matrix(~sex-1,data=D))}
  if ('income' %in% covariates) {D=cbind(D,model.matrix(~income-1,data=D))}
  if ('ethnicity' %in% covariates) {D=cbind(D,model.matrix(~ethnicity-1,data=D))}
  if ('education' %in% covariates) {D=cbind(D,model.matrix(~education-1,data=D))}
  D=cbind(D,model.matrix(~condition.w2-1,data=D))
  
  # intercept + covariates 
  controls=cbind(1,D$`topic.w2contractor pay`,D$`topic.w2farm payments`,D$`topic.w2foreign aid and abortion`,D$`topic.w2foreign workers`,D$`topic.w2gun violence research`,D$topic.w2LGBT,D$`topic.w2policing and military surplus`,D$`topic.w2public lands`,D$`topic.w2Russian sanctions`,D$`topic.w2student loans`,D$topic.w2trade,D$`topic.w2water rules`,D$`topic.w2weapon sales to Saudi Arabia`,D$topic.w2wildlife,D$president.w2Trump,D$condition.w2congress,D$condition.w2order)
  if ('topic.approve.w1' %in% covariates) {controls=cbind(controls,D$topic.approve.w1Disapprove,D$`topic.approve.w1Somewhat disapprove`,D$`topic.approve.w1Neither approve or disapprove`,D$`topic.approve.w1Somewhat approve`,D$`topic.approve.w1Approve`,D$`topic.approve.w1Strongly approve`)}
  if ('pres.approve.w1' %in% covariates) {controls=cbind(controls,D$pres.approve.w1Disapprove,D$`pres.approve.w1Somewhat disapprove`,D$`pres.approve.w1Neither approve or disapprove`,D$`pres.approve.w1Somewhat approve`,D$`pres.approve.w1Approve`,D$`pres.approve.w1Strongly approve`)}
  if ('manip.done.w1' %in% covariates) {controls=cbind(controls,D$manip.done.w1Disagree,D$`manip.done.w1Somewhat disagree`,D$`manip.done.w1Neither agree nor disagree`,D$`manip.done.w1Somewhat agree`,D$`manip.done.w1Agree`,D$`manip.done.w1Strongly agree`)}
  if ('trumpvote.w1' %in% covariates) {controls=cbind(controls,D$trumpvote.w1)}
  if ('sex' %in% covariates) {controls=cbind(controls,D$sexFemale)}
  if ('age' %in% covariates) {controls=cbind(controls,D$age)}
  if ('income' %in% covariates) {controls=cbind(controls,D$`income2: $25,000 to $50,000`,D$`income3: $50,001 to $75,000`,D$`income4: $75,001 to $100,000`,D$`income5: More than $100,001`)}
  if ('ethnicity' %in% covariates) {controls=cbind(controls,D$ethnicityBlack,D$`ethnicityNative American`,D$`ethnicityOther/Decline to State`,D$ethnicityWhite)}
  if ('education' %in% covariates) {controls=cbind(controls,D$`educationSome College or Vocational`,D$`educationB.A. or B.S.`,D$`educationPost-graduate or Higher`)}
  
  # new data to simulate over
  d.success.opp=cbind(controls,0,1,0,0,0)
  d.fail.opp=cbind(controls,0,1,1,0,1)
  
  d.success.cop=cbind(controls,1,0,0,0,0)
  d.fail.cop=cbind(controls,1,0,1,1,0)
  
  d.success.ind=cbind(controls,0,0,0,0,0)
  d.fail.ind=cbind(controls,0,0,1,0,0)
  
  # simulations
  p.win.opp=apply((1/(1+exp(-as.matrix(d.success.opp)%*% t(betas)))),1,as.vector) 
  p_win.opp=apply(p.win.opp,1,mean)
  p.fail.opp=apply((1/(1+exp(-as.matrix(d.fail.opp)%*% t(betas)))),1,as.vector) 
  p_fail.opp=apply(p.fail.opp,1,mean)
  
  p.win.cop=apply((1/(1+exp(-as.matrix(d.success.cop)%*% t(betas)))),1,as.vector) 
  p_win.cop=apply(p.win.cop,1,mean)
  p.fail.cop=apply((1/(1+exp(-as.matrix(d.fail.cop)%*% t(betas)))),1,as.vector) 
  p_fail.cop=apply(p.fail.cop,1,mean)
  
  p.win.ind=apply((1/(1+exp(-as.matrix(d.success.ind)%*% t(betas)))),1,as.vector) 
  p_win.ind=apply(p.win.ind,1,mean)
  p.fail.ind=apply((1/(1+exp(-as.matrix(d.fail.ind)%*% t(betas)))),1,as.vector) 
  p_fail.ind=apply(p.fail.ind,1,mean)
  
  # store
  results[1,1]<- quantile(p_fail.opp-p_win.opp,.0083)
  results[1,2]<- quantile(p_fail.opp-p_win.opp,.025)
  results[1,3]<- mean(p_fail.opp-p_win.opp)
  results[1,4]<- quantile(p_fail.opp-p_win.opp,.975)
  results[1,5]<- quantile(p_fail.opp-p_win.opp,.9916)
  results[1,6]<- 'Opposition'
  results[1,7]<- covariates[1]
  
  results[2,1]<- quantile(p_fail.cop-p_win.cop,.0083)
  results[2,2]<- quantile(p_fail.cop-p_win.cop,.025)
  results[2,3]<- mean(p_fail.cop-p_win.cop)
  results[2,4]<- quantile(p_fail.cop-p_win.cop,.975)
  results[2,5]<- quantile(p_fail.cop-p_win.cop,.9916)
  results[2,6]<- 'Copartisans'
  results[2,7]<- covariates[1]
  
  results[3,1]<- quantile(p_fail.ind-p_win.ind,.0083)
  results[3,2]<- quantile(p_fail.ind-p_win.ind,.025)
  results[3,3]<- mean(p_fail.ind-p_win.ind)
  results[3,4]<- quantile(p_fail.ind-p_win.ind,.975)
  results[3,5]<- quantile(p_fail.ind-p_win.ind,.9916)
  results[3,6]<- 'Independents'
  results[3,7]<- covariates[1]
  
  return(results)
}

other_outcomes=data.frame(lb.MC=numeric(12),lb.95=numeric(12),me=numeric(12),ub.95=numeric(12),ub.MC=numeric(12),grp=character(12),dv=character(12))
other_outcomes[1:3,1:7]=outcome_comparisons_party(t3a_c.bin,dt)
other_outcomes[4:6,1:7]=outcome_comparisons_party(p3a_c.bin,dt)
other_outcomes[7:9,1:7]=outcome_comparisons_party(v3a_c,dt)
other_outcomes[10:12,1:7]=outcome_comparisons_party(g3a_c.bin,dt)

# manip for figure
final_outcomes=rbind(other_outcomes,outcomes)

final_outcomes$lb.MC=as.numeric(final_outcomes$lb.MC)
final_outcomes$lb.95=as.numeric(final_outcomes$lb.95)
final_outcomes$me=as.numeric(final_outcomes$me)
final_outcomes$ub.95=as.numeric(final_outcomes$ub.95)
final_outcomes$ub.MC=as.numeric(final_outcomes$ub.MC)

final_outcomes$dv=as.factor(final_outcomes$dv)
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='trumpvote.w2'] <- 'Trump\nVote'
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='topic.approve.w2.bin'] <- 'Topic\nApproval'
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='pres.approve.w2.bin'] <- 'Job\nApproval'
levels(final_outcomes$dv)[levels(final_outcomes$dv)=='manip.done.w2.bin'] <- 'Gets things\nDone'

f2 = ggplot(final_outcomes[final_outcomes$dv != 'Gets things\nDone', ]) +
  geom_hline(yintercept = 0,
             colour = gray(1 / 2),
             lty = 2) +
  geom_linerange(
    aes(
      x = factor(
        dv,
        levels = c("Trump\nVote", "Job\nApproval", "Topic\nApproval")
      ),
      ymin = lb.MC,
      ymax = ub.MC
    ),
    lwd = 1.1,
    position = position_dodge(width = 1 / 2)
  ) +
  geom_pointrange(
    aes(
      x = factor(
        dv,
        levels = c("Trump\nVote", "Job\nApproval", "Topic\nApproval")
      ),
      y = me,
      ymin = lb.95,
      ymax = ub.95
    ),
    lwd = 2,
    size = 2.5,
    position = position_dodge(width = 1 / 2),
    shape = 21,
    fill = "WHITE"
  ) +
  geom_text(
    aes(
      x = factor(
        dv,
        levels = c("Trump\nVote", "Job\nApproval", "Topic\nApproval")
      ),
      y = me,
      label = round(me, 2)
    ),
    size = 3,
    position = position_dodge(width = 1 / 2)
  ) +
  facet_wrap( ~ grp, ncol = 4, scales = 'free_x') +
  coord_flip() + theme_bw() + xlab('') + ylab('') +
  theme(
    text = element_text(family = 'Palatino'),
    legend.position = 'bottom',
    legend.title = element_text(size = 9),
    strip.background = element_rect(
      color = "black",
      fill = "#ffffff",
      linewidth = 1.5,
      linetype = "blank"
    )
  ) +
  guides(color = guide_legend(
    override.aes = list(shape = 19, size = 1),
    reverse = T,
    title.position = "right"
  ))

ggsave("figures-tables/FIG-H2-1.pdf",
       plot = f2,
       height = 3,
       width = 7.2,
       units = 'in') 


f2_si = ggplot(final_outcomes[final_outcomes$dv == 'Gets things\nDone', ]) +
  geom_hline(yintercept = 0,
             colour = gray(1 / 2),
             lty = 2) +
  geom_linerange(aes(x = grp, ymin = lb.MC, ymax = ub.MC),
                 position = position_dodge(width = 1 / 2)) +
  geom_pointrange(
    aes(
      x = grp,
      y = me,
      ymin = lb.95,
      ymax = ub.95
    ),
    lwd = 2,
    size = 2.5,
    position = position_dodge(width = 1 / 2),
    shape = 21,
    fill = "WHITE"
  ) +
  geom_text(aes(x = grp, y = me, label = round(me, 2)),
            size = 3,
            position = position_dodge(width = 1 / 2)) +
  coord_flip() + theme_bw() + xlab('') + ylab('') +
  theme(
    text = element_text(family = 'Palatino'),
    legend.position = 'bottom',
    legend.title = element_text(size = 9),
    strip.background = element_rect(
      color = "black",
      fill = "#ffffff",
      linewidth = 1.5,
      linetype = "blank"
    )
  ) +
  guides(color = guide_legend(
    override.aes = list(shape = 19, size = 1),
    reverse = T,
    title.position = "right"
  ))

ggsave("figures-tables/FIG-H2-A-1.pdf",
       plot = f2_si,
       height = 3,
       width = 7.2,
       units = 'in') 


# Tabular outputs for appendix

stargazer(t3_c.bin,p3_c.bin,v3_c,g3_c.bin,digits=3,
          out = 'figures-tables/h2-tabular-output.txt',
          keep=c('condition','outcome'),
          omit.stat = c('ll','aic'),
          omit.table.layout = 'n',
          font.size='tiny',
          column.sep.width='3pt',
          column.labels=c('\\shortstack{Topic\\\\Approval}', '\\shortstack{Job\\\\ Approval}', '\\shortstack{Trump\\\\ Vote}', '\\shortstack{Gets things\\\\ Done}'),
          covariate.labels = c('Congress','Executive Order','Failure'),
          colnames = F, 
          dep.var.labels.include =F,
          model.numbers=F,
          title='{\\bf The public punishes presidents for failing to deliver.} Reports logistic regression coefficients and conventional standard errors with binary dependent variables indicating either approval, voting for incumbent, or positive agreement with question; simulated marginal effect estimates based on observed case approach from these models are reported in the leftmost panel of Figure \\ref{fig:h2}, and in Figure \\ref{fig:h2a}; all models condition on Wave 1 value of dependent variable, and also include topic and president factor variables, along with partisanship, age, income, sex, ethnicity, and education.',
          label='tab:h2',
          header=F)


stargazer(t3a_c.bin,p3a_c.bin,v3a_c,g3a_c.bin,digits=3,
          out = 'figures-tables/h2-tabular-output-2.txt',
          keep=c('condition','copartisanship','outcome'),
          omit.stat = c('ll','aic'),
          omit.table.layout = 'n',
          font.size='tiny',
          column.labels=c('\\shortstack{Topic\\\\Approval}', '\\shortstack{Job\\\\ Approval}', '\\shortstack{Trump\\\\ Vote}','\\shortstack{Gets things\\\\ Done}'),
          column.sep.width='3pt',
          colnames = F,
          covariate.labels=c('Congress','Exec. Order','Copartisan','Opposition','Failure','Copartisan X Failure','Opposition X Failure'),
          dep.var.labels.include =F,
          model.numbers=F,
          title='{\\bf Copartisans, the Opposition, and Independents punish presidents for failing to deliver.} Reports logistic regression coefficients and conventional standard errors with binary dependent variables indicating either approval, voting for incumbent, or positive agreement with question; simulated marginal effect estimates based on observed case approach from these models are reported in the right three panels of Figures \\ref{fig:h2} and \\ref{fig:h2a}; all models condition on Wave 1 value of dependent variable and also include topic and president factor variables, along with age, income, sex, ethnicity, and education.',
          label='tab:h2a',
          header=F)


stargazer(t3_c, p3_c, g3_c,
          digits = 3, 
          out = 'figures-tables/h2-tabular-output-ordered.txt',
          keep=c('condition','outcome'),
          omit.stat = c('ll','aic'),
          omit.table.layout = 'n',
          font.size='tiny',
          column.sep.width='3pt',
          covariate.labels = c('Congress','Executive Order','Failure'),
          column.labels=c('\\shortstack{Topic\\\\Approval}', '\\shortstack{Job\\\\ Approval}', '\\shortstack{Gets things\\\\ Done}'),
          colnames = F, 
          dep.var.labels.include =F,
          model.numbers=F,
          title='{\\bf The public punishes presidents for failing to deliver (Ordered Logits).} Reports ordered logistic regression coefficients and conventional standard errors with 7-point Likert dependent variables indicating ``Strongly disagree,`` ``Disagree,`` ``Somewhat disagree,`` ``Neither agree nor disagree,`` ``Somewhat agree,`` ``Agree,`` ``Strongly agree;`` all models condition on Wave 1 value of dependent variable, and also include topic and president factor variables, along with partisanship, age, income, sex, ethnicity, and education.',
          label='tab:h2_ordered',
          header = F)


stargazer(t3a_c, p3a_c, g3a_c,
          digits = 3, 
          out = 'figures-tables/h2-tabular-output-ordered-2.txt',
          keep=c('condition','copartisanship','outcome'),
          column.labels=c('\\shortstack{Topic\\\\Approval}', '\\shortstack{Job\\\\ Approval}','\\shortstack{Gets things\\\\ Done}'),
          omit.stat = c('ll','aic'),
          omit.table.layout = 'n',
          colnames = F,
          covariate.labels=c('Congress','Exec. Order','Copartisan','Opposition','Failure','Copartisan X Failure','Opposition X Failure'),
          dep.var.labels.include =F,
          model.numbers=F,
          font.size='tiny',
          column.sep.width='3pt',
          label='tab:h2a_ordered',
          header=F,
          title='{\\bf Copartisans, the Opposition, and Independents punish presidents for failing to deliver (Ordered Logits).} Reports ordered logistic regression coefficients and conventional standard errors with 7-point Likert dependent variables indicating ``Strongly disagree,`` ``Disagree,`` ``Somewhat disagree,`` ``Neither agree nor disagree,`` ``Somewhat agree,`` ``Agree,`` ``Strongly agree;`` all models condition on Wave 1 value of dependent variable and also include topic and president factor variables, along with age, income, sex, ethnicity, and education.')



